rm(list = ls())

library(sf)
library(ggplot2)
library(tigris)
library(haven)
library(maps)

#########################################

# set file paths (pull from Stata file 00_main.do)
raw_path          <- Sys.getenv("RAW")
intermediate_path <- Sys.getenv("INTERMEDIATE")
for_analysis_path <- Sys.getenv("FOR_ANALYSIS")
temp_path         <- Sys.getenv("TEMP")
output_path       <- Sys.getenv("OUTPUT")

geo_data <- file.path(raw_path, "geo_data")
ipums_puma_2000 <- file.path(geo_data, "ipums_puma_2000")

# read in data on SNAP retailer accessibility by PUMA
food_access <- read_dta(file.path(intermediate_path, "puma_food_access_R.dta"))
food_access <- food_access[, c("STATEFIP", "PUMA", "LA1and10")]

# read in IPUMS shape file from https://usa.ipums.org/usa/volii/boundaries.shtml
shape <- st_read(file.path(ipums_puma_2000, "ipums_puma_2000.shp"))
states <- map_data("state")

# need to omit Alaska and Hawaii
merged_access <- merge(shape, food_access, by = c("STATEFIP", "PUMA"), all.x = T)
merged_access <- merged_access[!(merged_access$STATEFIP == 15),]
merged_access <- merged_access[!(merged_access$STATEFIP == "02"),]

# generate and save figure
access_map <- ggplot() +
  geom_sf(data = merged_access, aes(fill = LA1and10)) +
  scale_fill_continuous(name = "Low Access Probability",
                        low = "lightblue", high = "darkblue", 
                        na.value = "grey75") +
  theme_void()

ggplot2::ggsave(filename = file.path(output_path, "access_map.png"), plot = access_map, device = "png", dpi = 1200)
